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In  this  effort  we  will  design  and  develop  a  multi  spectral  sensing  imaging  and  communi¬ 
cations  system  based  on  a  compressive  sensing  (CS)  architecture  that  can  computationally 
adapt  dynamically,  in  real  time,  to  a  set  of  optimal  spectral  band(s)  needed  to  address  the 
desired  tactical  missions.  The  effort  integrates  the  development  and  implementation  of  state 
of  the  art  signal  processing  algorithms  and  communications,  especially  compressive  sensing 
algorithms,  with  coded  aperture  multispectral  sensors  and  optical  design  concepts,  to  attain 
a  compact  spectrally  agile  sensor  system.  The  system  under  development  is  composed  of 
three  subsystems: 

(1)  Relay/dispersion  and  focal  plane  array  (FPA)  compressed  measurements. 

(2)  Imaging  and  coded  apertures. 

(3)  Analog  source /channel  coding  subsystem.  In  this  quarter,  the  following  progress  has 
been  attained  in  each  of  these  thrusts: 

For  the  period  of  performance  ending  June  2010  the  following  has  been  accomplished  under 
this  contract  in  each  of  the  subsystems: 

1  Relay  Compressed  Measurements 

1.1  Installed  a  monochromator-based  setup  to  acquire  calibration 
cubes  needed  by  the  image  reconstruction  algorithm 

The  goal  of  the  proposed  imager  is  to  obtain  information  about  the  imaging  scene  by  mea¬ 
suring  the  intensity  of  the  light  at  each  spatial  location  at  different  wavelengths  which  helps 
to  obtain  detailed  information  of  the  object.  Thus,  in  this  part  of  the  effort,  we  installed  an 
elaborate  setup  involving  a  monochromator  along  with  other  optical  parts  to  acquire  calibra¬ 
tion  cubes  which  could  be  used  in  the  reconstruction  of  the  image.  The  image  reconstruction 
algorithm  needs  to  be  trained  before  being  used  in  the  image  reconstruction  experiment.  The 
training  process  mainly  involves  the  acquisition  of  a  calibration  cube,  which  is  the  a  priori 
knowledge  of  the  spatially  shifted  information  of  the  photomask  in  different  spectral  chan¬ 
nels.  The  calibration  cube  can  be  different  when  different  lenses/prisms  are  used  to  build  the 
system.  To  experimentally  realize  such  a  training  process  (acquiring  the  calibration  cube), 
we  built  a  spectrally  adjustable  illumination  source  using  a  monochromator  and  a  Xenon 
(Xe)  lamp.  The  monochromator  used  in  our  experiments  is  supplied  by  the  Oriel  Corpora¬ 
tion  (Oriel  model  number:  77200).  Using  suitable  gratings  and  splits,  this  monochromator 
can  achieve  a  spectral  resolution  of  1  nm.  The  output  spectrum  of  the  monochromator  can 
be  adjusted  by  rotating  the  blazing-grating  in  the  monochromator.  Thus,  various  spectrum 
ranges  can  be  obtained  with  different  gratings.  For  our  monochromator,  the  spectrum-range 
is  440-670nm  (visible).  A  Xe  lamp  (Newport  model  number:  66477)  is  used  as  the  input  to 
the  monochromator.  Figure  1  shows  the  monochromator  setup. 
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Figure  1:  Monochromator  setup.  The  output  spectrum  of  the  monochromator  can  be  ad¬ 
justed  by  rotating  blazing-grating.  The  spectrum-range  is  440-670nm  (visible). 

1.2  Acquisition  of  the  calibration  cube 

The  key  components  in  the  process  of  obtaining  the  calibration  cube  are  the  monochromator, 
photomask,  prism,  relay  lens  and  a  CCD  camera.  Particularly,  the  monochromator  serves 
as  the  source;  photomask  is  used  for  the  random  aperture  code  while  the  prism  is  used  to 
introduce  spectral  dispersion.  A  relay  lens  relays  the  image  from  the  plane  of  the  coded 
aperture  to  the  CCD  which  captures  the  dispersed  image.  Moreover,  we  used  prisms  of 
different  apex  angles  to  acquire  the  calibration  cube,  including  an  equilateral  prism  (Apex 
angle:  60),  a  right  angle  prisms  (Apex  angle:  45),  a  wedge  prism  (Apex  angle:  18),  and 
the  double-amici  prism  mentioned  in  our  last  quarterly  report  in  order  to  analyze  their 
performances.  Figure  2  shows  the  experimental  setup  used  to  acquire  the  calibration  cube. 

Figure  3  shows  an  example  of  a  calibration  cube  obtained  using  the  equilateral  prism. 
The  spectral  resolution  of  the  monochromator  is  lnm.  Therefore,  for  the  given  spectral 
range,  there  are  231  slices  in  the  calibration  cube.  We  show  21  of  them  in  Fig.  3.  The 
adjacent  spectral  images  in  this  figure  have  a  10  nm  difference  in  monochromator-emission 
wavelength. 

In  addition  to  equilateral  prisms,  we  also  used  prisms  with  apex  angles  of  450  and  180  to 
acquire  the  calibration  cube.  In  particular,  Fig.  4  compares  the  calibration  cubes  obtained 
with  different  prisms.  Generally  speaking,  the  addition  of  prism  introduces  aberrations  to 
the  image  formation  process.  For  instance,  we  can  see  the  edges  of  the  calibration  cubes 
have  some  curvatures.  The  aberrations  are  caused  by  the  geometric  shape  of  the  prism.  The 
prism  has  different  thicknesses  at  different  heights  with  respect  to  the  apex  angle.  When  the 
light  propagates  through  the  prism  at  different  height/angle,  its  propagation  will  be  affected 
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Figure  2:  Experimental  setup  used  to  acquire  the  calibration  cube  using  an  equilateral  prism 
with  apex  angle  (apex  angle:  60) 
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Figure  3:  Calibration  cube  obtained  using  an  equilateral  prism  (apex  angle:  60). 

by  different  thicknesses  of  the  glass  material.  Since  the  aberrations  are  considered  in  the 
training  process,  the  image  reconstruction  algorithm  can  minimize  the  quality  degradation 
caused  by  the  aberrations.  In  Fig.  4,  we  can  see  that  calibration  cubes  acquired  using  prisms 
of  smaller  apex  angles  have  apparently  less  held  curvature. 

2  Code  Aperture  Agile  Spectral  Imaging  System  (CAASI) 

Recently,  the  Coded  Aperture  Snapshot  Spectral  Imaging  (CASSI)  architecture  has  made 
it  possible  to  implement  CS  in  spectral  imaging.  CASSI  is  indeed  a  remarkable  imaging 
architecture  that  has  been  studied  extensively  in  [1,  3,  4],  The  single-shot  CASSI  architec¬ 
ture,  however,  may  use  excessive  compression  to  represent  spectrally  rich  image  cubes  under 
surveillance,  leading  in  some  cases  to  low  quality  image  reconstructions  as  well  as  low  spec¬ 
tral  resolution.  The  coding  and  reconstruction  algorithms  in  CASSI  are  also  rigid  because 
the  entire  spectral  image  cube  is  reconstructed  at  once.  The  new  CAASI  aims  to  overcome 
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Figure  4:  Calibration  cubes  obtained  using  the  prisms  of  different  apex  angles. 
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Figure  5:  Representation  of  the  CAASI  parts.  The  £x  code  aperture  in  the  traditional  CASSI 
system  is  replaced  by  a  spatial  light  modulator. 


this  difficulties  and  provides  the  flexibility  to  recover  a  specific  subset  of  spectral  bands  with 
controllable  SNR  or  reconstruction  time. 

The  compressive  code  aperture  single  shot  spectral  imaging  system  is  depicted  in  Fig. 
5  [1].  The  coding  is  realized  by  the  coded  aperture  T(x,y)  as  applied  to  the  image  source 
density  /0(x,  y\  A)  where  (x,  y )  are  the  spatial  coordinates  and  A  is  the  wavelength  resulting  in 
the  coded  field  f\  (x,  y,  A).  The  coded  density  is  spectrally  dispersed  by  a  dispersive  element 
before  it  impinges  on  the  focal  plane  array  (FPA)  as  f-2(x,y,  A), 

/2(x,  y,  X)  =  j  J  T(x',  y’)f0(x',  y',  A )h(x’  -  aX  -  x,  y’  -  y)dx'dy’  (1) 

where  T{x' ,  y ')  is  the  transmission  function  representing  the  code  aperture,  h(x'— 0(A)—  x,  y' — 
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y)  is  the  optical  impulse  response  of  the  system,  and  a\  is  the  dispersion  induced  by  the  prism 
assuming  a  linear  dispersion.  The  compressive  measurements  across  the  FPA  are  realized  by 
the  integration  of  the  held  /2(x,  y,  A)  over  the  detector’s  spectral  range  sensitivity.  The  source 
density  can  be  written  in  discrete  form  as  (F k)nm  where  n,m  index  the  spatial  coordinates, 
and  k  determines  the  kth  spectral  plane.  Following  the  mathematical  model  in  [4],  the  coding 
is  realized  by  an  aperture  pattern  (T)nm.  The  compressed  sensing  measurements  at  the  focal 
plane  array  can  be  written  in  the  following  matrix  form: 

L 

(G)  nm  ^  ^(^k^n,m+k(^^n,m+k  +  un  (2) 

k=  1 


where  unm  represents  white  noise,  L  is  the  number  of  spectral  bands  and  where  n  G 
{1,2,  and  m  G  {1,2,...,M}  index  the  pixels  on  the  detector.  The  expression  in 

(2)  can  be  expressed  as 

g  =  Hf  +  uj  =  HTd  +  u  (3) 

where  g  and  f  are  a  vector  representation  of  G  and  F  respectively.  H  is  the  projection  matrix, 
T  is  a  Kronecker  basis  representation,  and  9  is  a  sparse  coefficient  vector  representing  f.  If 
the  aperture  code  pattern  is  fixed  and  only  one  snap-shot  is  detected,  the  resultant  spectral 
imager  is  the  so-called  single  disperser,  or  single-shot,  CASSI  architecture  [2],  In  this  case, 
the  entire  3-D  multispectral  image  cube  is  compressed  into  a  single  2-D  compressive  image 
measurement  at  the  FPA.  The  compression  factor  rj,  defined  as 


V 


M  +  L-l 
M  x  L 


(4) 


indicates  the  grade  of  compressive  measurements.  As  rj  — »  1  the  linear  set  of  equations  in  (3) 
became  severely  underdetermined  causing  the  CASSI  measurements  to  be  highly  compressed, 
and  in  consequence  the  accurate  reconstruction  of  the  original  information  is  difficult.  On  the 
contrary,  as  r)  — >  0  little  compression  is  attained  making  it  easier  to  recover  the  information 
embedded  in  G.  The  spectral  image  cube  f  can  be  reconstructed  by  solving  the  optimization 
problem  f  =  \k{argmin0,  ||G  —  H\k#/|||  +  t||0,||i}  where  r  >  0  is  a  regularization  parameter 
that  balances  the  conflicting  tasks  of  minimizing  the  least  square  of  the  residuals  while,  at 
the  same  time,  yielding  a  sparse  solution[5].  In  our  work,  we  have  obtained  an  equivalent 
model  equation  of  CASSI  system  in  alternative  form  as 


gq  =  DCqfq,  (5) 

where  D  represents  the  operation  of  the  dispersive  element,  and  Cq  is  the  matrix  represen¬ 
tation  of  the  qth  row  of  the  code  aperture.  The  vector  representation  in  (5)  characterizing 
CASSI,  will  be  very  useful  in  our  project. 


3  Analog  Source/channel  Coding  Subsystem 

In  this  work,  we  use  analog  joint  source-channel  coding  for  the  transmission  of  data  samples 
at  high  rates.  Here,  we  study  the  performance  of  such  a  system  in  comparison  with  optimized 
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Figure  6:  Block  diagram  of  an  N:1  bandwidth  compression  analog  joint  source-channel  coded 
system. 

capacity  approaching  digital  Bit  Interleaved  Coded  Modulation  (BICM)  schemes.  The  aim 
is  to  determine  under  which  circumstances  it  is  more  convenient  to  use  each  communication 
setup.  We  measure  the  performance  of  both  analog  and  digital  systems  in  terms  of  Signal- 
to-Distortion  Ratio  (SDR)  versus  Channel  Signal-to-Noise  Ratio  (CSNR)  when  Gaussian 
and  Laplacian  distributed  sources  are  transmitted.  Notice  that  Gaussian  and  Laplacian 
sources  are  of  interest  for  reasons  beyond  their  theoretical  significance.  In  fact,  images 
represented  in  transform  domains  (e.g.,  wavelets)  can  be  accurately  modeled  using  these 
distributions,  which  makes  the  analog  joint  source-channel  coding  system  very  relevant  for 
image  communications  [12],  We  show  that  analog  transmission  performs  better  than  the 
digital  scheme  with  a  much  lower  encoding  and  decoding  complexity. 

3.1  Analog  joint  source-channel  coded  system 

The  distortion  between  source  symbols  X  =  and  decoded  symbols  X  =  is 

N 

calculated  according  to  the  MSE,  defined  as  MSE  =  A  Yh  I \xi  ~  xi\  |2-  Consequently,  the 

i=  1 

system  performance  can  be  measured  in  terms  of  the  output  SDR  with  respect  to  the  CSNR, 
with  SDR  defined  as  SDR  =  10  log  (1/MSE),  where  the  source  symbols  are  normalized  to 
unit  mean  power.  Given  N  and  K,  the  theoretival  limit  (OPTA)  is  calculated  by  equating 
the  rate  distortion  function  to  the  AWGN  channel  capacity  [8].  For  example,  for  Gaussian 
sources  we  have  IV log  (1  /MSE)  =  /Flog  (1  +  1/ofJ. 

Figure  6  shows  the  block  diagram  of  an  N:1  analog  joint  source-channel  coding  system 
where  the  source  generates  blocks  of  B  i.i.d.  symbols  that  are  encoded  into  B/N  channel 
symbols.  Without  loss  of  generality,  we  assume  a  source  distribution  with  zero  mean  and 
unit  variance  and  also  that  the  mean  transmit  power  is  equal  to  one.  In  this  paper,  we 
focus  on  memoryless  Gaussian  and  Laplacian  sources.  A  particular  type  of  parameterized 
space-filling  continuous  curves,  called  spiral-like  curves  [8]- [10],  can  be  used  to  encode  the 
X  =  (xl,x2)  source  samples.  For  the  case  of  2  :  1  compression  ( i.e.,N  =  2),  they  are 
formally  defined  as 

I  xe,i  =  signed)  £ sin(9 ) 

\xe,2  =  ^cosiO)  for  6»  <G  77 

where  A  is  the  distance  between  two  neighboring  spiral  arms,  and  6  is  the  angle  from  the 
origin  to  the  point  Xg  =  (xitV,X2,v)  on  the  curve.  Therefore,  each  pair  of  source  samples,  xl 
and  x2,  represent  a  specific  point  in  77  2  that  is  matched  to  the  closest  point  Xg  =  (xi^,x2,w) 
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Figure  7:  Block  diagram  of  a  Bit  Interleaved  Coded  Modulation(BICM)  system  over  AWGN 
channel. 

on  the  spiral.  The  angle  from  the  origin  to  that  point  on  the  spiral,  6  ,  will  be  the  channel 
symbol  for  xl  and  x2,  i.e. 

6  =  Ma  (A")  =  ar  g  min  {(xi  ±  (A/V)  6>sin0)2  +  (x2  —  (A/7r)  OcosO)2  }  (7) 

Since  our  goal  is  the  minimization  of  the  MSE,  the  bidimensional  space  has  to  be  filled  by 
the  spiral  in  the  best  possible  way  for  every  CSNR  value.  On  one  hand,  by  changing  the  A 
value,  we  manage  to  optimize  this  matching  and  to  improve  the  system  performance.  On 
the  other  hand,  it  is  possible  to  achieve  higher  compression  rates  (i.e.  N:l)  by  extending  (1) 
to  generate  more  complex  curves  [14],  [15]. 

The  next  step  consists  in  defining  an  invertible  function  of  0-with  the  corresponding 
normalization  factor  to  ensure  the  transmit  power  constraint.  In  [8,  10,  16]  the  invertible 
function  Ta{6)  =  9a  with  a  =  2,  was  proposed.  However,  as  shown  in  [11],  the  system 
performance  can  be  improved  if  a  and  A  are  numerically  optimized  for  each  different  CSNR 
value.  Therefore,  the  channel  symbol  is  Ta  (jfj  / ,  where  y'y  is  the  normalization  factor. 

In  summary,  the  received  symbol  y  at  the  decoder  can  be  expressed  as  y  =  Ta  (Ma  (A"))  + 
nyj 7,  where  n  is  Gaussian  noise.  Given  a  received  symbol  y,  MMSE  decoding  is  performed  at 
the  receiver  to  calculate  an  estimation  of  the  corresponding  source  symbol.  Optimal  MMSE 
decoding  can  be  expressed  as 

A mmse  =  E  {X\y}  =  I X p  (X \y)  dX  =  (1/p  (y))  J  X p  (y\X)  p  (A)  dX  (8) 

where  the  mapping  function  Ma  (•)  is  used  to  obtain  the  conditional  probability  p(y\ X).  Note 
that  the  integral  in  (8)  can  only  be  calculated  numerically  because  Ma  (•)  is  discontinuous 
and  highly  non-linear.  To  do  so,  X  is  first  discretized  using  a  uniform  step  and  a  mapped 
value  is  calculated  for  each  discretized  point  according  to  (7).  As  a  result,  we  obtain  a 
discretized  version  of  p(y\X).  Next,  p( X)  is  also  computed  for  each  point,  and  thus  the 
calculation  of  the  integral  is  reduced  to  multiplicative  and  additive  operations.  Since  this 
discretization  does  not  depend  on  the  received  symbol,  it  is  calculated  once  off-line  and 
stored  in  the  decoder.  Although  in  this  paper  we  focus  on  2:1  systems,  the  proposed  system 
can  be  readily  modified  to  adapt  the  compression  rate  from  N:1  to  N:K  [17]. 

3.2  Bit  interleaved  coded  modulation  (BICM) 

Figure  7  shows  the  block  diagram  of  a  digital  BICM  system.  We  assume  a  discrete-time 
source  that  produces  Gaussian  and  Laplacian  independent  and  identically  distributed  (i.i.d.) 


real  valued  analog  symbols  Xa.  These  continuous  samples  are  mapped  to  a  discrete  set 
of  values  using  an  optimum  Q-level  scalar  quantizer.  Both,  the  quantization  levels  and 
the  partition  regions  can  be  obtained  using  the  well-known  Lloyd-Max  algorithm  [21,  22], 
Although  better  performance  could  be  obtained  with  vector  quantization,  we  discarded  this 
possibility  to  keep  the  overall  quantizing  complexity  at  a  low  level.  Next,  the  quantized 
discrete-time  symbols  are  converted  into  a  binary  representation  using  a  suitable  source 
encoder.  Again,  among  the  many  existing  source  encoding  methods,  we  decided  to  use 
Huffman  encoding  because  it  is  a  simple  algorithm  and  approaches  the  source  entropy.  The 
input  alphabet  to  the  Huffman  encoder  is  made  up  of  the  Q-levels  of  the  scalar  quantizer, 
i.e.,  no  grouping  is  performed  prior  to  the  encoding.  The  length  of  the  average  Huffman 
codeword  used  to  represent  each  source  sample  will  be  denoted  by  Lm.  The  output  bit 
sequence  is  encoded  with  a  rate  r  capacity  approaching  channel  encoder.  Due  to  their  low 
encoding  and  decoding  complexity,  we  decided  to  use  Irregular  Repeat  Accumulate  (IRA) 
codes  [23].  Finally,  the  channel  encoded  bits  are  modulated  using  a  real- valued  M-PAM 
constellation.  Notice  that  an  interleaver  between  the  channel  encoder  and  the  modulator  is 
not  strictly  necessary  since  we  are  transmitting  over  an  AWGN  channel.  A  one-dimensional 
PAM  constellation  has  been  chosen  to  keep  the  same  signaling  as  the  analog  system  in  Section 
II  where  real  analog  encoded  symbols  are  transmitted.  The  constellation  size  M  limits  the 
maximum  attainable  data  rate  over  the  channel.  We  chose  the  value  M  =  256  which  allows 
the  transmission  of  a  maximum  of  8  bits  per  channel  use,  a  transmission  rate  high  enough 
for  the  comparison  carried  out  in  the  ensuing  section.  Higher  values  of  M  could  be  used  but 
they  yield  to  extremely  complex  PAM  constellations  that  are  not  feasible  in  practice.  At 
the  receiver,  an  optimum  detector  calculates  the  LLRs  of  the  transmitted  bits  and  passes 
them  to  the  sum-product  IRA  decoding  algorithm.  After  a  maximum  number  of  decoding 
iterations,  the  resulting  bits  are  hard-decoded  and  dequantized  to  the  corresponding  levels. 
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